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Recently, the prospects for amorphous phases of 
graphene (a-G) have been explored computationally. 
Initial models were flat, and contained odd-member 
rings, while maintaining three-fold coordination and sp 2 
bonding. Upon relaxation, puckering occurs, and may be 
traced to the existence of pentagons, in analogy with the 
situation for fullerenes. In this work, we systematically 
explore the inherent structures with energy close to the 
flat starting structure. As expected, the planar symmetry 
can be broken in various ways, which we characterize for 
800-atom model of a-G, always using local basis 



density functional techniques. The classical normal 
modes of various structural models are discussed, with 
an emphasis on imaginary modes indicating the evolu- 
tion from flat to puckered. We also discuss very low en- 
ergy conformational fluctuations akin to those seen pre- 
viously in amorphous silicon [8], and reflect on the na- 
ture of the amorphous "ground state" within a network 
of fixed topology. For completeness, high energy modes 
were also computed, and are found to be associated with 
strained parts of the network. 
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1 Introduction A key feature of matter in a disor- 
dered state is the existence of large numbers of conforma- 
tions with essentially degenerate energies, which may also 
be accessible with small energy cost. This is in contrast 
with crystals, which possess long range order, few poly- 
morphs and typically a deep energy minimum and large en- 
ergy barrier. The ground state is thus sharply defined, and 
the only low energy excitations are phonons - transitions to 
other structures are prohibited. Where realistic models of 
disordered systems are concerned, few attempts have been 
made to quantitatively characterize the number, energet- 
ics and proximity (in the sense of barrier) of these states. 
In his inherent structure formulation of statistical mechan- 
ics, Stillinger argued that the number of minima scales like 
N\exp(aN)[l], where N is the number of atoms in the 
model, and a is a positive system-dependent constant, a 
was estimated to be around 0.8 in a monatomic liquid[2], 
and flexible organic molecules exhibit larger a, as in the 
fragile glass former ortho-terphenyl where a is around 



13. 14 [3]. For temperatures well below the melting point, 
these local energy minima are denoted inherent structures. 
Each inherent structure is confined to its own basin, which 
refers to a set of configurations can be mapped to the in- 
herent structure by potential energy minimization!!]. 

To further motivate this work, consider the following 
gedanken experiment. Consider a sequence of molecules 
with TV atoms. It is well known that as N increases, the 
number of minima accessible to the molecule also increase, 
and such conformations are extensively studied in chem- 
istry^] [5] [6] [7]. While it is unlikely that a rigorous theory 
quantifying these minima as a function of N can be formu- 
lated, it is clear from computer experiments that the num- 
ber of minima grow drastically with N. For most molecu- 
lar systems it is difficult to be certain that the global mini- 
mum structure has been found in a simulation as there are 
so many metastable minima in which the system can be 
trapped. In this paper, we are concerned with the even more 
intractable problem of characterizing the minima, or poten- 
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tial energy surface of a disordered condensed matter sys- 
tem in both two-dimensional (2D) and three-dimensional 
(3D) cases. For a 3D crystal, if one introduces small ran- 
dom distortions and relaxes the disturbed system, it returns 
to exactly the same structure. This is not true for disor- 
dered systems. For 3D, we find there exits a continuum of 
metastable minima for a-Si, in which a number of tiny dis- 
tortions of bond angles (and to a much lesser extend bond 
lengths) yield a distinct energy-degenerate conformation, 
which reveals the existence of an extraordinarily flat poten- 
tial energy landscape (PEL). Fedders and Drabold showed 
that starting from a well-relaxed a-Si:H model, quenching 
a set of "snapshots" during a constant-T MD simulation 
never returns to the exact initial state, instead they fall into 
minima that are topologically equivalent (eg. with the same 
network connectivity), but with small variations in bond 
angles and bond lengths[8]. We find that the behavior in a- 
Si is consistent with a-Si:H. The case of a-G is different 
in the following sense. Like a-Si, we find a continuum of 
essentially energy degenerate states in which minute (but 
"real") variations in bond angles and bond lengths are dis- 
played. Of course, such states retain identical connectivity 
since their energies are identical to within a few micro eV. 
However, a-G also exhibits a variety of local energy min- 
ima associated with different puckering. These structures 
usually have similar energies (within ~ lOmeV), but a 
significant barrier separating them. Thus, the picture that 
emerges of the a-G energy landscape is a variety of inher- 
ent structures (with varying puckering) with slightly vary- 
ing energies but substantial barriers between them, and in 
each of these basins is associated with a particular puck- 
ered state, an ambiguously defined minimum with small 
variations in bond angles and bond lengths accessible as 
we describe in detail below in Sec. 4.2. 

Beside the work described above on ultra-low energy 
excitations, we also discuss other phenomena peculiar to 
a-G. There has been intensive study in understanding the 
properties of crystalline graphene, but little is securely un- 
derstood about amorphous phases. Recent electron bom- 
bardment experiments have revealed the existence of amor- 
phous graphene[17] [18]. Clear images of regions of amor- 
phous graphene have been published by Meyer[26]. In pre- 
vious work, we observed that planar amorphous graphene 
is extremely sensitive to out-of-plane distortions [9] . Simi- 
lar behavior has been verified in amorphous graphene by 
experiment and other calculations [11] [12] and also in 
silicon nanosheets [13]. We have found different initial 
conditions in breaking the planar symmetry of amorphous 
graphene model lead to distinct puckered states after re- 
laxation [9]. These states exhibit little difference in topo- 
logical properties, ie. ring statistics and coordination num- 
ber. However, the total energy differences between these 
metastable states are around 0.02eV per atom, and the full 
width of the puckering along the original normal direction 
is around 6 — 7 A. 




a) Top view of 800-atom crystalline graphene. 




b) Top view of 800-atom amorphous graphene due to He and Thorpe 
laxed by us. 

Figure 1 Comparison between crystalline and amorphous phases 
of graphene. 

A natural complement to these studies is an exploration 
of low-frequency classical normal modes. These modes 
turn out to be rather delocalized. As presented for two level 
systems, the tunneling between two equilibrium states trig- 
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gers a number of low-energy excitations [15] [16]. Details 
are discussed in Sec. 4.3. In Sec. 5, we summarize our find- 
ings about the degenerate states and localized imaginary-, 
low- and high-frequency vibrational modes of amorphous 
graphene. 

2 Models To approach this problem computation- 
ally, we employ an 800-atom amorphous graphene model 
(800 a-g) due to He and Thorpe generated by introducing 
Stone- Wales defects into a perfect honeycomb lattice and 
a WWW annealing scheme[27]. This model has perfect 
threefold coordination with varying concentration of 5, 6 
and 7 member rings [10] and is a practical realization of the 
continuous random network (CRN) concept proposed by 
Zachariason[28]. The comparison between crystalline and 
amorphous graphene is shown in Fig. 1. We have relaxed 
the 800-atom a-g model with an accurate ab-initio code, 
while small rearrangements occurred and planar symmetry 
was preserved. We found that by very slightly breaking 
the planar symmetry (by randomly moving each atom by 
~ 0.01 A) and performing a structural optimization, the 
resulting minima were always puckered[9] . 

The amorphous Si model we employ is a realistic 
64-atom model (64 a-Si), generated by Barkema and 
Mousseau using a modified form of WWW algorithm [29]. 
This model has perfect fourfold coordination and within 
the limitations of its small size, to our knowledge is not in 
significant contradiction to any experiment. 

3 Procedure Our calculations are performed with an 
ab-initio program SIESTA[30], using pseudopotentials and 
the Perdew-Zunger parameterization of the local-density 
approximation (LDA) with a single-^ basis and Harris- 
Functional at a constant volume. For earlier simulations 
the details of method are described in [19] and [20]. To 
investigate the nature of minima on the potential energy 
surface, we employ a method proposed by Fedders and 
Drabold[8] similar to the conformational space annealing 
approach mentioned in [21], which has been used in lo- 
cating and predicting low-energy conformations of var- 
ious proteins[22][23][24][25]. First, starting with a per- 
fectly relaxed model (in our case 800 a-g and 64 a-Si), 
we run a sequence of four parallel simulations. We let 
the network evolve for 8. Ops at four different mean tem- 
peratures of 20K, 500K, 600K or 900K. The target tem- 
peratures are achieved by velocity rescaling. From these 
simulations,we drew instantaneous snapshots of the con- 
figuration by saving the instantaneous configuration ev- 
ery 0.1bps, then rapidly quench these to OK and let the 
quenched configurations evolve around OK for a few ps to 
find the metastable minimum (or inherent structure) asso- 
ciated with the initial snapshots. 

To investigate the topological changes between these 
minima (quenched configurations from snapshots), we use 
two autocorrelation functions as defined by Fedders and 



Drabold[8]: 

A0(t u t 2 ) = J2m(h)-0 z (t 2 )) 2 /N^ (1) 

i 

and 

Ar{h,t 2 ) = E(( r ^i) -n(t2)) 2 /N) k 2 (2) 

In Eq. 1 the index i runs through all the bond angles where 
Oi is the i th bond angle. In Eq. 2 the index i runs over all 
nearest-neighbor pairs where is the i th is the distance 
between a pair. All the atoms in quenched supercells are 
threefold, making the definitions of Oi and well defined. 
The times t\ and t 2 refer to the quenched snapshots. These 
autocorrelation functions provide a close view of how ther- 
mal MD simulations induce transitions between different 
energy basins. 

4 Discussion We break the discussion into three 
parts. First, the nature of pentagonal puckering, confor- 
mation fluctuations and finally an analysis of the classical 
vibrational modes. 



i 1 1 1 1 1 1 1 1 1 1 1 r 




Figure 2 Color. Correlation between the total energy per atom 
and magnitude of puckering for constant temperature MD simu- 
lations. The zero total energy refers to the total energy of original 
flat 800 a-g model. 

4.1 Symmetry Breaking As we have shown in [9], 
the original flat 800-atom a-g model is exceedingly sensi- 
tive to transverse distortion, and then loses planar symme- 
try lowering the total energy of the supercell. In every case, 
even at T = 20K, the planar symmetry breaks and the sys- 
tem puckers: thermal disorder is sufficient to induce puck- 
ering. Fig. 2 shows the relation between the total energy 
of the system and maximum separation of atoms along the 
normal direction (magnitude of puckering) in constant-T 
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MD simulations at the indicated temperatures. In the lan- 
guage of PEL, starting from the flat 800 a-g, four MD sim- 
ulations overcome tiny energy barriers and take a down-hill 
path to regions with lower energy. Thus the flat 800 a-g 
model can be considered as an exceptionally shallow basin 
on the PEL. The barrier to puckering from flat is a few mi- 
cro eV for this Hamiltonian. 



0.5 



A 0.25 




1 1 1 1 

— Quenched minima from 500K MD 

— ■ Quenched minima from 600K MD 
. • — Quenched minima from 900K MD 




Figure 3 Color. Time variation of two autocorrelation functions. 
This figure shows autocorrelation functions of A0(ti,t2) and 
100Z\r(£i, £2) for £1 = 6. Ops and £2 varying from 6.0 to 7.95ps. 
The temperatures are 500K, 600K and 900K. The functions ap- 
pear to be continuous. 
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Figure 4 Color. Total energy distribution functions of quenched 
supercells from MD runs under 20K, 500K, 600K and 900K. The 
total energy of original flat 800 a-g is considered as OeV. Dis- 
tinct structures correspond to different puckered states, broad- 
ening within each major peak from conformational variations. 
Three major peaks from MD runs at 600K, 500K and 900K are 
labeled as 1, 2 and 3 respectively. 



4.2 Conformational Fluctuations The quenching 
procedures at sequential timesteps yield basins on the PEL. 
Here we show the calculations at average temperatures of 
500K, 600K and 900K. These results of MD runs at dif- 
ferent temperatures exhibit consistency with each other. 
Fig. 3 shows how these two autocorrelation functions vary 
with time of snapshots. Since temperatures of all the MD 
simulations achieve equilibrium after 6.0 ps, here the au- 
tocorrelation functions are calculated with t\ = 6. Ops. It 
appears A6(ti, £2) and Ar(ti,t2) from three MD runs at 
different temperatures are qualitatively similar. They in- 
crease linearly and approach a constant. The continuity of 
the curves in Fig. 3 suggests that there is a continuum of 
states, accessible albeit structurally varying only in very 
modest ways. 




Figure 5 Color. Side view of two quenched configurations. Gray 
balls and sticks show the configuration from 900K MD, and blue 
lines represent the one from 500K. 

The total energy distributions of all the quenched su- 
percells from MD runs at different temperatures are shown 
in Fig. 4. For temperatures other than 20K, the total energy 
distributions exhibit several peaks. The minor peaks in Fig. 
4 correspond with the annealing process of MD runs. The 
major peaks (labeled 3, 2 and 1 in Fig. 4) are derived from 
different puckering configurations sampled in the process 
of equilibration to constant T. Correspondingly as shown 
in Fig. 3, the fluctuations of autocorrelation functions (af- 
ter thermal equilibrium is reached) reach an asymptotic 
state after 6. Ops. Each of the three peaks in Fig. 4 demon- 
strates a basin on the PEL of a-G. These nearly degenerate 
quenched equilibrium states are trapped in distinct basins 
on the PEL, and quenched minima within one basin form a 
continuum metastable state around inherent structures. De- 
tails of variations in bond angles, bond lengths and total en- 
ergies of these metastable states are shown in Table 1 . The 
total energy variation between the basins corresponding to 
the major three peaks are averaged as 1.405 x 10 _3 eV, 
which is one order of magnitude higher than the energy of 
fluctuations within a basin. In spite of their different energy 
scales, these quenched configurations belonged to distinct 
basins share identical local bonding. The only difference is 
that they pucker in distinct ways, as shown in Fig. 5. 

For comparison, we repeat parallel calculations using 
64 a-Si model quenched from MD runs at 20K, 300K, 
500K and 600K. The variations of autocorrelation func- 
tions are shown in Fig. 6. The results are in agreement with 
Fedders and Drabold[8]. We see for a-Si systems, there ex- 
ists one general basin on the PEL (for a particular network 
connectivity), and the paths lowering the total energy on 
the PEL will eventually go into this basin, leading to in- 
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Table 1 Average value and standard deviation of E to t /N a tom, Ar(t\ , £2) and A0(ti , £2) of quenched configurations from MD runs in 
the time period from 3.6 to 8.0 ps, where £1 = 1.05ps. 



T(K) 


Etot /N a tom (eV) 


O 'Etot /N atom 


(eV) 


Ar(l,t 2 )(A) 


(JAr (A) 


A0(l,t 2 )(o) 


OAQ (0) 


500K 


-9.595 x 10" 2 


8.855 x 10" 


-5 


4.907 x 10-2 


2.551 x 10" 4 


3.736 


1.814 x 10" 2 


600K 


-9.435 x 1(T 2 


9.917 x 10" 


-5 


7.503 x 10~ 2 


1.973 x 10~ 4 


4.181 


6.464 x 10" 2 


900K 


-9.716 x 10" 2 


2.091 x 10" 


-4 


2.516 x 10" 2 


8.08 x 10" 3 


1.725 


0.146 




time (P s ) frequency (cm" 1 ) 



Figure 6 Color. Time variation of two autocorrelation func- 
tions for a-Si. This figure shows autocorrelation functions of 
AQ{t\,t<2) and 100Z\r(£i, £2) for £1 = 6.6ps and £2 varying 
from 6.6 to 8.0 ps. The temperatures are 20K, 300K, 500K and 
600K. The results are similar to [8]. 



herent structures with minor changes in bond angles and 
lengths, and analogous energy scale. 

Comparison between results of a-G and a-Si suggests 
that PEL of 3D system (a-Si) is smooth and inherent struc- 
tures are contained in one general basin. a-G is similar 
within one puckered state. 




Figure 7 Color. Side view of pucker-up and -down 800 a-g mod- 
els. Gray balls and sticks illustrate pucker-up model, and pucker- 
down supercell is represented by blue lines. 

For a-G, an MD run at higher temperature (annealing) 
can overcome the energy barrier between basins and reach 
an inherent structure with lower energy. Also by investigat- 
ing into the correlation between topology and energy scale 
of these quenched supercells, it is revealed that lower total 
energy (stabler state) is related with more puckered con- 
figurations with small variation in bond lengths and angles 
from the original flat 800-atom model. 



Figure 8 Color. Vibrational density of states (VDOS) of 800 
crystalline graphene, pucker-up and -down a-g models. Note the 
distractive feature at uj = 1375cm -1 for a-G. 



4.3 Classical Normal Modes To investigate the vi- 
brational modes in these supercells, we perform calcula- 
tions of dynamical matrix eigenvalues and eigenvectors of 
original flat 800 a-g model, and two quenched configura- 
tions with certain region puckering along opposite direc- 
tion, designated "pucker-up" and "pucker-down" models, 
as shown in Fig. 7. The dynamical matrix was constructed 
by moving each atom along six directions by 0.04 Bohr. We 
also perform the phonon calculation for a 800-atom crys- 
talline graphene model. The vibrational density of states 
(VDOS) of 800 crystalline model, pucker-up and -down 
800 a-g models are shown in Fig. 8. The VDOS result of 
crystalline graphene shows good agreement with published 
calculation[31]. In Fig. 8 at a frequency near 1375 cm -1 , 
the spectrum of crystalline graphene has a minimum. In 
contrast the spectrum of two puckered supercells achieve 
a local maximum. Thus Raman scattering experiments are 
expected to provide a way to distinguish crystalline and 
amorphous graphene. There is no difference in the spec- 
trum between pucker-up and -down 800 a-g models. 

In the original flat 800-atom a-g model, the eigenvec- 
tors with imaginary eigenvalues have large components 
along the normal direction of the graphene plane (at least 
four orders of magnitude higher than x and y components). 
These imaginary-frequency modes are localized on pen- 
tagons in the network: two example are shown in Fig. 9. 
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flat 800 a-g Freq = -283.843 cm" 1 flat 800 a "§ Frec l = "272.661 cm 1 




a) Imaginary-frequency mode in flat 800 a-g model with fre- b) Imaginary-frequency mode in flat 800 a-g model with fre- 
quency uj — 283.843icra -1 (referred as im-model). quency uj — 272.661icm~ 1 (referred as im-mode2). 




c) Side view of the region where im-model is originally lo- d) Side view of the region where im-mode2 is originally lo- 
calized, calized. 



Figure 9 Color. Two examples of imaginary-frequency modes in flat 800 a-g model. The contour plot represents the component of 
eigenvector along the direction transverse to the plane. 



30 



— ■ pucker-up 800 a-g model 
■ — pucker-down 800 a-g model 



s 
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Figure 10 Color. Temperature dependence of C(T) of pucker-up 
and -down 800 a-g models. 



As shown in Fig. 9c and 9d, these imaginary-frequency 
modes are localized near structures that lead to pucker- 
ing. As shown in [9], it is pentagons that lead to puckering 
and symmetry breaking. Thus these imaginary-frequency 



modes are an indicator of the instability of the flat 800 a-g 
model. 

In the puckered models, we observe modes with a low 
frequency, around 14-20 cm -1 , reminiscent of "floppy 
modes" proposed by Phillips and Thorpe[34][35]. The 
structures of these low-frequency modes are quite com- 
plex, as shown in Fig. 11. These modes are rather ex- 
tended, and have significant weight on pentagonal puck- 
ered regions and large rings, analogous to what Fedders 
and Drabold have seen in a-Si:H[8]. The observed energy 
scale of these low-frequency modes is around a few meV, 
almost half of the lowest frequency of the acoustic phonon 
modes in a crystalline graphene with same size. As pro- 
posed in the theory of "two-level systems", there exists 
a vast distribution of low-energy excitations, caused by 
tunneling of atoms between nearly degenerate equilibrium 
states[15][16]. Goldstein pointed out that the dynamics 
could be separated into two categories: vibrational mo- 
tion about a minimum on PEL and transitions between 
minima[36]. Then these low-frequency modes might be 
triggered by transitions between degenerate minima within 
one basin on the PEL. As shown in Sec. 4.2, the energy 
variations between minima within one basin is in the order 
of 10 _4 eV, and the energy difference between basins is 
in the order of 10 _3 eV. The energy scales of these low- 
frequency modes (~ meV) are sufficient to drive confor- 
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10 20 30 40 



a) Low-frequency mode in pucker-down 800 ce-g model with fre- 
quency uj — 19.167cm -1 . 




b) Low frequency mode in pucker-up 800 a-g model with fre- 
quency uj — 17.318cm -1 . 

Figure 11 Color. Examples of low-frequency modes in pucker- 
down and -up 800 a-g models. The contour plots represent the 
intensity of eigenvectors on each atom. The blue atoms illustrate 
the "puckering-most" atoms, and the green atoms represent "flat" 
atoms. 



mational fluctuations, but not high enough to overcome the 
energy barrier between different basins (quenched states) 
on the PEL. 

In the high-frequency domain, there also exist highly 
localized high-frequency modes in the puckered configura- 
tions. These modes are triggered by the pentagonal defects 
and are highly localized, as shown in Fig. 12. This result is 
consistent with Biswas et al. [37] and Fedders et al. [38], 
who have shown strain and topological defects are active 
at highest frequencies. 



Pucker- Down Freq = 1642.073 cm 
Intensity plot 




a) High-frequency mode in pucker-down 800 a-g model with fre- 
quency uj — 1642.073cm -1 . 

Pucker-Up Freq = 1629.252 cm" 1 
Intensity plot 




10 20 30 40 

x 



b) High-frequency mode in pucker-up 800 a-g model with fre- 
quency uj = 1629.252cm -1 . 

Figure 12 Color. Examples of localized high-frequency modes 
in pucker-down and -up 800 ce-g models. The contour plots rep- 
resent the intensity of eigenvectors on each atom. The blue atoms 
illustrate the "puckering-most" atoms, and the green atoms repre- 
sent "flat" atoms. 



We also compute the specific heat C (T) using VDOS 
information [3 9]: 

rE max F p E/k B T 

c{T) = 3R l <^> V/^-i)' 9(£)<i£ 

(3) 

where g(E) is normalized VDOS. For room temperature 
(300K), specific heat of flat, pucker-up and -down 800 a- 
g models are 25.151, 18.879 and 18.702 JK^mol- 1 re- 
spectively. The temperature dependence of C(T) is shown 
in Fig. 10. This is presumably an academic result as it is 
currently hard to imagine an experiment for C (T) for this 
2D system. 
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5 Conclusion In conclusion, we have found that a- 
G has a rich and interesting energy landscape. We observe 
distinct energy scale of basins (~ lOmeV) associated with 
different puckered configurations and then within such a 
configuration, an ambiguous energy minimum with a con- 
tinuum of bond angles and bond lengths with energy scale 
(~ few/ieV) and a nearly flat PEL. Within a given puck- 
ered configuration, this continuum is much like what was 
seen for a-Si:H in 1996[8]. 

Vibrational calculations reveal the existence of local- 
ized imaginary-frequency modes in flat 800 a-g model. 
These modes are localized on pentagons and play the key 
role in losing planar symmetry and forming pentagonal 
puckering structures. We find delocalized low-frequency 
phonon modes, similar to floppy modes, which have sub- 
stantial weight on defects and share the same energy scale 
as the energy difference between adjacent basins on the 
PEL. Thus these low-frequency modes are triggered by the 
transition between adjacent energy minima. Some high- 
frequency modes are detected and highly localized on 
puckered regions and large rings. 

Acknowledgements Here we want to strongly acknowl- 
edge Dr. M. F. Thorpe at Arizona State University and his former 
student Harry He, who have prepared the 800-atom a-G model. 
This work is supported by NSF Grant No. DMR 09-03225. 

References 

[1] F. H. Stillinger, Science 267, 1935 (1995). 
[2] D. C. Wallace, Phys. Rev. E 56, 4179 (1997). 
[3] F. H. Stillinger, J. Phys. Chem. B 102, 2807 (1998). 
[4] F. H. Stillinger, and T A. Weber, Phys. Rev. A 28, 2408 
(1983). 

[5] F. H. Stillinger, Phys. Rev. E 59, 48 (1999). 
[6] R. Schilling, Physica D 216, 157 (2006). 
[7] J. P. K. Doye, and D. J. Wales, J. Chem. Phys. 116, 3777 
(2002). 

[8] P. A. Fedders, and D. A. Drabold, Phys. Rev. B 53, 3841 
(1996). 

[9] Y. Li, F. Inam, A. Kumar, M. Thorpe, and D. A. Drabold, 

Phys. Stat. Sol. B 248, 2082 (2011). 
[10] V. Kapko, D. A. Drabold, and M. F. Thorpe, Phys. Status 

Solidi B 247, 1197-1200 (2010). 
[11] M. T Lusk, and L. D. Carr, Phys. Rev. Lett. 100, 175503 

(2008). 

[12] E. Holmstrom, J. Fransson, O. Eriksson, R. Lizarraga, B. 
Sanyal, S. Bhandary, and M. Katsnelson, Phys. Rev. B 84, 
205414(2011). 

[13] E. C. Anota, A. B. Hernandez, M. Castro, and G. H. Coco- 
letzi, J. Mol. Model. 18, 2147-2152 (2012). 

[14] W A. Phillips, Amorphous Solids, Low Temperature Prop- 
erties, Springer, Berline 1981. 

[15] P. W Anderson, B. I. Halperin, and C. M. Varma, Philos. 
Mag. 25, 1 (1972). 

[16] W A. Phillips, J. Low Temp. Phys. 7, 351 (1972). 

[17] J. Kotakoski, J. C. Meyer, S. Kurasch, D. Santos-Cottin, U. 
Kaiser, and A. V. Krasheninnikov, Phys. Rev. B 83, 245420 
(2011). 



[18] J. Kotakoski, A. V. Krasheninnikov, U. Kauser, and J. C 

Meyer, Phys. Rev. Lett. 106, 105505 (2011). 
[19] D. A. Drabold, R. Wang, S. Klemm, O. F. Sankey, and J. D. 

Dow, Phys. Rev. B 43, 5132 (1991). 
[20] D. A. Drabold, P. A. Fedders, and M. P. Grumbach, Phys. 

Rev. B 54, 5480 (1996). 
[21] D. J. Wales, Energy Landscapes with Applications to 

clusters, Biomolecules and Glasses, Cambridge University 

Press, New York, 2003. 
[22] J. Lee, H. A. Scheraga, and S. Rackovsky, J. Comput. Chem. 

18,1222(1997). 
[23] J. Lee, H. A. Scheraga, and S. Rackovsky, Biopolymers 46, 

103 (1998). 

[24] J. Lee, and H. A. Scheraga, Int. J. Quant. Chem. 75, 255 
(1999). 

[25] A. Liwo, J. Lee, D. R. Ropoll, J. Pillardy, and H. A. Scher- 
aga, Proc. Natl. Acad. Sci. USA 96, 5482 (1999). 

[26] J. C. Meyer, C. Kisielowski, R. Erni, M. D. Rossell, M. F. 
Crommie, and A. Zettl, Nano. Lett. 8, 3582 (2008). 

[27] F. Wooten, K. Winer, and D. Weaire, Phys. Rev. Lett. 54, 
1392 (1985). 

[28] W H. Zachariasen, J. Am. CHem. Soc. 54, 3841 (1932). 
[29] P. Biswas, R. Atta-Fynn, and D. A. Drabold, Phys. Rev. B 

69, 195207 (2004). 
[30] E. Artacho, E. Anglada, O. Dieguez, J. D. Gale, A. Garcia, 

J. Junquera, R. M. Martin, P. Ordejon, J. M. Pruneda, D. 

Sanchez-Protal, and J. M. Soler, J. Phys. Condens. Matter 

20, 064208 (2008). 
[31] F. Liu, P. Ming, and J. Li, Phys. Rev. B 76, 064120 (2007). 
[32] A. Hashimoto, K. Suenaga, A. Gloter, K. Urita, and S. Li- 

jima, Nature (London) 430, 870 (2004). 
[33] J. Kotakoski, A. V. Krasheninnikov, U. Kaiser, and J. C. 

Meyer, Phys. Rev. Lett. 106, 105505 (2011). 
[34] J. C. Phillips, J. Non-Cryst. Solids 43, 37 (1981). 
[35] M. F. Thorpe, J. Non-Cryst. Solids 57, 355 (1983). 
[36] M. Goldstein, J. Chem. Phys. 51, 3728 (1969). 
[37] R. Biswas, M. Bouchard, W A. Kamitakahara, G. S. Grest, 

and C. M. Soukoulis, Phys. Rev. Lett. 60, 2280 (1988). 
[38] P. A. Fedders, D. A. Drabold, and S. Klemm, Phys. Rev. B 

45, 4048 (1992). 
[39] S. M. Nakhmanson, and D. A Drabold, Phys. Rev. B 61, 

5376 (2000). 



Copyright line will be provided by the publisher 



